clear all;
%% Code explanation


%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

L_W1=64/50.6;  % sample 1: 3Sn/8PbTe
L_W2=48.5/52.9; % sample 2: 3Sn/6PbTe


ratio=0.5;
%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20180122 3Sn_6PbTe 3Sn_8PbTe outofplanedata\'); %

filenames={'20180122_0uW_wroots';'20180122_20uW_wroots';'20180122_50uW_wroots';...
           '20180122_100uW_wroots';'20180122_150uW_wroots';'20180122_200uW_wroots';...
           '20180123_50uW_onlyrot';'20180122_250uW_wroots';'20180123_150uW_onlyrot';...
           '20180122_0uW_onlyrot';'20180123_200uW_onlyrot';'20180123_300uW_onlyrot';...
           '20180123_450uW_onlyrot';'20180123_600uW_onlyrot';'20180123_750uW_onlyrot';'20180124_900uW_onlyrot';...
           '20180124_1050uW_onlyrot'};


Boffset=0.011;
NT=length(filenames);
colorset=jet(NT);
figure

Tmc=zeros(NT,1);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT,2);
Rc=zeros(NT,2);
for kk=1:NT
   
    file1=[root,filenames{kk},'_perp_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f','commentStyle','#');
        data1=data1';
        fclose(fid);
 
          B=data1{1}-Boffset;
          R1=(data1{2}/L_W1+0.009)/10;
          R2=(data1{3}/L_W2+0.02)/10;
          T=T_RuOx(data1{4}*100);
          Tmc(kk,1)=mean(T);
          Tmc_err(kk,1)=3*std(T);
          
          subplot(1,2,1),plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
          hold on
          subplot(1,2,2),plot(B,R2,'Color',colorset(kk,:)','LineWidth',1);
          hold on

          
          
          indexB=find(B>1.2);
          p1=polyfit(B(indexB),R1(indexB),1);
          p2=polyfit(B(indexB),R2(indexB),1);
          
          index1=find((R1-ratio*(p1(1)*B+p1(2)))<=0);
          if isempty(index1)==0
              Hc2(kk,1)=B(index1(1));
              Rc(kk,1)=R1(index1(1));
          else
              Hc2(kk,1)=NaN;
              Rc(kk,1)=NaN;
          end
          
          index2=find((R2-ratio*(p2(1)*B+p2(2)))<=0);
          if isempty(index2)==0
              Hc2(kk,2)=B(index2(1));
              Rc(kk,2)=R2(index2(1));
          else
              Hc2(kk,2)=NaN;
              Rc(kk,2)=NaN;
          end
          title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[B R2 T Tmc(kk,1)*ones(length(B),1)];
       
         xlswrite('FigS4A.xlsx',AA,kk);
       xlswrite('FigS4A.xlsx',title,kk);    
 
end
index1=find(isnan(Hc2(:,1))==0);
index2=find(isnan(Hc2(:,2))==0);
p11=polyfit(Hc2(index1,1),Rc(index1,1),1);
p22=polyfit(Hc2(index2,2),Rc(index2,2),1);
subplot(1,2,1),plot(B,p11(1)*B+p11(2),'-k');
              hold on
subplot(1,2,2),plot(B,p22(1)*B+p22(2),'-k');
              hold on
Hc21=Hc2(index1,1);
Hc22=Hc2(index2,2);
Tmc1=Tmc(index1);
Tmc2=Tmc(index2);
N1=length(Tmc1);
M1=N1-2;
N2=length(Tmc2);
M2=N2-2;
pp1=polyfit(Tmc1(M1:N1),Hc21(M1:N1),1);
pp2=polyfit(Tmc2(M2:N2),Hc22(M2:N2),1);

Tmc10=abs(pp1(2)/pp1(1));
Tmc20=abs(pp2(2)/pp2(1));

xi2=sqrt(6.626E-34/4/pi/1.6E-19/(abs(pp2(2))))*1E9;

Tp=(0:0.002:Tmc20)';

figure
plot([Tmc1' Tmc10]',[Hc21' 0]','dr',[Tmc2' Tmc20]',[Hc22' 0]','sb',Tp,pp2(1)*Tp+pp2(2),'-k');
hold on





